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Nf = 2 + 1 dynamical Wilson quark simulation toward the physical point 



1. Introduction 



The first lattice QCD calculation of hadron masses in 1981 revealed its potential ability to non- 
perturbatively evaluate the physical quantities in the strong interaction from first principles. Since 
then, the history of lattice QCD has been a succession of enduring efforts to control the major 
systematic errors due to finite lattice size, finite lattice spacing, quenching and chiral extrapolation. 
Thanks to recent progress of simulation algorithms and increasing availability of computational 
resources, we are about to bring all the above systematic errors under control. This will allow us 
to establish whether or not QCD is the fundamental theory of the strong interaction by investi- 
gating the hadron spectrum, and further proceed to elucidate the fundamental issues of the strong 
interactions and the Standard Model. 

The previous CP-PACS/JLQCD project[|I|, ||] aimed at a full removal of the quenching effects 
by performing Nf = 2 + 1 lattice QCD simulations with the nonperturbatively 0(a)-improved Wil- 
son quark actional and the Iwasaki gauge actional on a (2 fm) 3 lattice at three lattice spacings. 
While we have succeeded in incorporating the dynamical strange quark effects by the Polyno- 
mial Hybrid Monte Carlo (PHMC) algorithm the lightest up-down quark mass reached with 
the HMC algorithm was 64 MeV corresponding to mjz/mp « 0.6, which required a long chiral 
extrapolation to the physical point at m n /m p w 0.18. 

The PACS-CS (Parallel Array Computer System for Computational Science) project^, [7], ||] is 
the successor to the CP-PACS/JLQCD project, which takes up the task that the latter has left off, 
namely simulation at the physical point to remove the ambiguity of chiral extrapolation. It employs 
the same quark and gauge actions as the CP-PACS/JLQCD project, but uses the PACS-CS computer 
with a total peak speed of 14.3 TFLOPS developed and installed at University of Tsukuba on 1 
July 2006. The up-down quark masses are reduced by employing the domain-decomposed HMC 
(DDHMC) algorithm with the replay trick proposed by Luscher[||, So far, we have reached 
the up-down quark mass of 6 MeV which yields the pion mass of about 210 MeV. We also improve 



the simulation of the strange quark part with the UV-filtered PHMC (UV-PHMC) algorithm! 1 1 



In this report we present simulation details and preliminary results which include the chiral 
analysis on the pseudoscalar meson masses and the decay constants with chiral perturbation theory, 
the light hadron spectrum and the p-%% mixing effects. Some algorithmic issues are also discussed. 
Selected topics on the light hadron spectrum and the ChPT analysis on the pseudoscalar meson 
masses and the decay constants are also reported in Refs.[12, l3fj. 



2. Simulation details 



2.1 Simulation parameters 

We employ the f?(a , )-improved Wilson quark action with a nonperturbative improvement co- 
efficient csw = 1.715|^] and the Iwasaki gauge action at j8 = 1.90 on a 32 3 x 64 lattice which is 
enlarged from 20 3 x 40 in the CP-PACS/JLQCD project to investigate the baryon masses. Simula- 
tion parameters are summarized in Table |l| We choose six combinations of the hopping parameters 
("lid) K s) based on the previous CP-PACS/JLQCD results. Among them the heaviest combination 
(K- ud) K- s ) = (0.13700,0.13640) in this work is the lightest one in the previous CP-PACS/JLQCD 
simulations, which enable us to make a direct comparison of the two results with different lattice 
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sizes. As for the strange quark, the hopping parameter K s = 0.13640 corresponds to the physi- 
cal point JC S = 0.136412(50) as estimated in the CP-PACS/JLQCD work[jl|, §]. This is the reason 
why all our simulations are carried out with K s = 0.13640, the one exception being the run at 
K s = 0.13660 and K" lK j = 0.13754 to investigate the strange quark mass dependence. 

In order to simulate the degenerate up-down quarks we employ the DDHMC algorithm, whose 
effectiveness for small quark mass region has already been shown in the Nf = 2 case[^, [l4|, [fj]. 
The characteristic feature of this algorithm is a geometric separation of the up-down quark de- 
terminant into the UV and the IR parts as a preconditioner of HMC, which is implemented by 
domain-decomposing the full lattices with small blocks. We choose 8 4 for the block size being less 
than (1 fm) 4 in physical units and small enough to reside within a computing node of the PACS-CS 
computer. There are two prominent points in the DDHMC algorithm. Firstly, communication be- 
tween the computing nodes is not required in calculating the UV part, which is a preferable feature 
for alleviating the problem of a widening gap between the processor floating point performance and 
the network communication bandwidth with parallel computers. Secondly, we can incorporate the 



multiple time scale integration scheme Q16Q to reduce the simulation cost efficiently. The relative 



magnitudes of the force terms are found to be 

H^gl|:||*uv||:||fiR||» 16:4:1, (2.1) 

where we adopt the convention ||M|| 2 = — 2tr(M 2 ) for the norm of an element M of the SU(3) 
Lie algebra. F g denotes the gauge part and Fjv.ir are for the UV and the IR parts of the up-down 
quarks. The associated step sizes for the forces are controlled by three integers A^o,i,2 : = 
z/N N\N 2 , <5tuv = x/N\Nz, <5t ir = f/A^ with % the trajectory length. The integers #0,1,2 are 
chosen such that 

5Tg||F g || SttjvII^uvII ~ <5tir||Fir||. (2.2) 



Taking account of the relative magnitudes of the forces in eq.(2.1 ) we find a larger value is allowed 
for 5tir compared to 5t g and 5?uv, which means that we need to calculate Fir less frequently 
in the molecular dynamics trajectories. Since the calculation of Fir contains the quark matrix 
inversion on the full lattice, which is the most time consuming part, this integration scheme saves 
the simulation cost remarkably. The values for #0,1.2 are listed in Table [jj where Nq and Ni are fixed 
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38.2(17.3) 
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(4,4,14) 
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20.9(10.2) 
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(4,4,20) 
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19.2(8.6) 




0.13660 


0.5 


(4,4,28) 
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10.3(2.9) 


0.13770 


0.13640 


0.25 


(4,4,16) 


180 


2000 


38.4(25.2) 


0.13781 


0.13640 


0.25 


(4,4,48) 


180 


350 


9.1(6.1) 



Table 1: Simulation parameters. MD time is the number of trajectories multiplied by the trajectory length 
T. Tim [P] denotes the integrated autocorrelation time for the plaquette. 
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Figure 1: Plaquette history (left) and normalized autocorrelation function (right) for (jc„d, ff s ) = 
(0.13727,0.13640). Horizontal lines in the left denote the average value of the plaquette with an error 
band. 



at 4 for all the hopping parameters, while the value of Af 2 is adjusted taking account of acceptance 
rate and simulation stability. 

For the UV-PHMC algorithm for the strange quark, the domain-decomposition is not imple- 
mented. Since we have found | |F S | | ~ | |^ir| |> the step size is chosen as 5t s = 5tir. The polynomial 
order for UV-PHMC, which is denoted by A^ po i y in Table [l], is adjusted to yield high acceptance rate 
for the global Metropolis test at the end of each trajectory. 

The inversion of the Wilson-Dirac operator D on the full lattice is carried out by the SAP 
(Schwarz alternating procedure) preconditioned GCR solver, where the preconditioning can be ac- 
celerated with the single-precision arithmetic whereas the GCR solver is implemented with the 



double precision[17]. We employ the stopping condition \Dx — b\/\b\ < 10 for the force calcula- 
tion and 10~ 14 for the Hamiltonian, which guarantees the reversibility of the molecular dynamics 
trajectories to high precision: \AU\ < 10~ 12 for the link variables and \AH\ < 10~ 8 for the Hamil- 
tonian at (fc ud ,K- s ) = (0.13781,0.13640). 

2.2 Plaquette history and autocorrelation time 

In Fig. [I] we show the plaquette history and the normalized autocorrelation function p(r) 
at (JCudjKs) = (0.13727,0.13640) as a representative case. The integrated autocorrelation time is 
estimated Ti n t [P] = 20.9(10.2) following the definition in Ref. [g] 



Tint(T) = ~ + 



I P(T) 



(2.3) 



0<T<W 



where the summation window W is set to the first time lag z that p(f) becomes consistent with 
zero within the error bar. In this case we find W = 119.5. The choice of W is not critical for 
estimate of Ti nt in spite of the long tail observed in Fig. [lj Extending the summation window, we 
find that Tj nt [P] saturates at Ti nt [P] ~ 25 beyond T = 200, which is within the error bar of the original 
estimate. For other hopping parameters we have found similar behaviors for the plaquette history 
and the normalized autocorrelation function. We hardly observe the quark mass dependence for 
Tnt [P] listed in Table [j] The statistics may not be sufficiently large to derive a definite conclusion, 
however. 
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Figure 2: Cost estimate of Nf = 2 + 1 QCD simulations with the HMC (solid line) and the DDHMC (red 
circle) algorithms at a = 0.1 fm with L = 3 fm for 100 independent configurations. Vertical dotted line 
denotes the physical point. 



3. Algorithmic issues 



3.1 Efficiency of DDHMC algorithm 

In order to discuss the efficiency of the DDHMC algorithm, it is instructive to compare with 
that of the HMC algorithm. We first recall an empirical cost formula for Nf = 2 QCD simulations 



with the HMC algorithm based on the CP-PACS results []18|]: 



cost [Tfl ops • years] =C 



#conf 
1000 



0.6 



m K /m p 



i 5 



3 fm 



0.1 fm 



n 7 



with C « 2.8. A strong quark mass dependence is found in the above formula: 1 / (m n /m p ) e behaves 
as in the leading term for the small quark mass region. This quark mass dependence is 

owing to the following three factors: The number of the quark matrix inversion is governed by 
the condition number which should be proportional to 1 /m u &; to keep the acceptance rate fixed we 
should take 8x oc m uc j for the step size in the molecular dynamics trajectories; The autocorrelation 



time of the HMC evolution shows 1 /m U( \ dependences in the CP-PACS run[|T9|]. 

To estimate the computational cost for Nf = 2 + 1 QCD simulations with the HMC algorithm, 
we assume that the strange quark contribution is given by half of eq.(3T) at m n /m p = 0.67 which 



is a phenomenologically estimated ratio of the strange pseudoscalar meson "m^" and 



nir 



2nv\ - 



mt 



0.67. 



(3-1) 



Since the strange quark is relatively heavy, its computational cost occupies only a small fraction 
as the up-down quark masses decrease. In Fig. ^| we draw the cost formula for the Nf = 2 + 1 
case as a function of m K jmp, where we take #conf=100, a=0.1 fm and L = 3 fm in eq.([3~l|) as a 
representative case. We observe a steep increase of the computational cost below m n /m p ~ 0.5. At 
the physical point the expected cost is 0(100) Tflops-years. 
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(N ,N h N 2 ) 


Npoly 


trajs. 


MD time 




Tint[#mult] 


0.5 


(4,4,6) 


130 


6000 


3000 


23.7(9.2) 


92(56) 


0.5/3 


(4,4,2) 


130 


18000 


3000 


18.6(5.9) 


42(21) 



Table 2: Parameters for T-dependence study. #mult denotes the number of multiplications of the Wilson- 
Dirac quark matrix on the full lattice. 

Now let us turn to the case of the DDHMC algorithm. The red symbol denotes the measured 
cost at fc u d=0.13781, 0.13770 with ?c s = 0.13640, which are the lightest two points in our simula- 
tion. The DDHMC algorithm show a remarkable improvement reducing the cost by 30 — 50 times 
in magnitude. The majority of this reduction arises from the multiple time scale integration scheme 
and the GCR solver accelerated by the SAP preconditioning with the single-precision arithmetic. 
Roughly speaking, the improvement factor is O(10) for the former and 3 — 4 for the latter. Note 
that the quark mass dependence is also tamed: Since we find that Ti nt [P] is independent of the 
quark masses, the cost is proportional to l/m^ d . Our results show a feasibility of simulations at the 
physical point with the 0(10) Tflops computer which is available at present. 

3.2 T dependence of DDHMC algorithm 

In the DDHMC algorithm a subset of of all link variables, which are referred to as the active 
link variables, are updated during the molecular dynamics evolution, while keeping other field 
variables fixed[]9]]. The fraction of the active link variables depends on the block size we choose. In 
our case of 8 4 it is only 37%. To ensure that all the link variables on the lattice should be updated 
equally on average we implement random gauge field translations at the end of every molecular 
dynamics trajectories following Ref. [g]. 

Our concern is that the DDHMC algorithm might have a long autocorrelation time due to 
the existence of fixed link variables during the molecular dynamics evolution. A possible way 
to reduce the effects of the fixed link variables is more frequent random gauge field translations. 
This is easily realized by making z shorter with 5t fixed. We investigate the T dependence of the 
DDHMC algorithm employing a smaller lattice size of 16 3 x 32 at (fc ud , K s ) = (0.13700,0.13640). 
Other parameters are summarized in Table |2[ 




MD time MD time 

Figure 3: Normalized autocorrelation functions for the plaquette (left) and the #mult (right) at (fc u d, K s ) = 
(0.13700,0.13640). Black (red) symbols denote the T = 0.5 (t = 0.5/3) case. 
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In Fig. |3] we show the normalized autocorrelation functions for the plaquette and the number 
of multiplications of the Wilson-Dirac quark matrix on the full lattice during a molecular dynamics 
trajectory. Black symbols denote the X = 0.5 case, while red ones are for the z = 0.5/3 case. We 
observe that the normalized autocorrelation functions for z = 0.5 have longer tails than z = 0.5/3 
before becoming consistent with zero: This is quantitatively checked by the integrated autocorrela- 
tion times Tint[P] and Tj nt [#mult] in Table ||: The z = 0.5/3 case has shorter autocorrelation times. 
Since the computational cost is the same for the z = 0.5 and the z = 0.5/3 cases in terms of MD 
time, we can conclude that the z = 0.5/3 case shows a better efficiency than the z = 0.5 case. 
Based on this study, albeit conducted at a relatively heavy quark mass (m K /m p « 0.6), we employ 
z = 0.25 for the production run at (fc ud , ?c s ) = (0.13781,0.13640) and (0.13770,0.13640), which 
is half of the trajectory length at other hopping parameters. 

3.3 Simulation stability 



In Refs. [14, 15] simulation stability was discussed based on the spectral gap distribution of 



the Wilson-Dirac operator for two-flavor lattice QCD simulations. The spectral gap is defined as 

pL = min{|A| | A is an eigenvalue of Q}, (3.2) 

where Q is the hermitian Wilson-Dirac operator Q = JsDw with Dw = {^/2){jj 1 (V* jl + V^,) — 
aV* V^} + mo- Important indices to characterize the distribution are its median p. and width a. 
The latter is defined as (v — u)/2, where [u,v] is the smallest range of /I which contains more than 
68.3% of the data. This is to avoid potentially large statistical uncertainties which might occur 
when data are not sufficiently sampled. Their chief findings are two points: The first one is that 
the median jl shows a good linear dependence on the current up-down quark mass m^ 1 and the 
magnitude of the slope is well described by Za empirically. The second one is that the width a 
scales as 



a— ~ 1 (3.3) 
a 

with V the four-dimensional volume in physical units. They also observe that the width a is roughly 
independent of the quark mass for the unimproved Wilson quark action, while it shows a trend to 
decrease with the mass for the improved one. 

This study was also applied to the Nf = 2 + 1 case in Ref. [Q] where we reported on our 
preliminary run on a 16 3 x 32 lattice preparing for the PACS-CS project. We observed jl oc m^ 1 
and found 0.5<a(\/V/a)<0.76 for 15 MeV < m^ 1 < 64 MeV, where o diminishes as the up- 
down quark mass decreases. 

The existence of a gap in the spectrum of the Wilson-Dirac operator allows us to simulate the 
light quarks efficiently. The authors in Ref. Jl4] ] propose a stability condition requiring jl > 3a to 
assure the existence of the gap. Let us apply this condition to our case. Assuming a(^/V /a) = 1 
we estimate a = 2.26 MeV using a = 0.09 fm and V = (2.8 fm) 4 which will be obtained later. 
By using the empirical relation jl ~ Zawj^* 1 we find Z J 4m^ ( J VI >6.8 MeV for the stability condition, 
which is heavier than the physical point. On the other hand, we found o(VV /a) < 1 in Ref. [S], 
indicating that the actual bound will be lower. Our runs toward the physical point should shed light 
on the actual bound of stability for our lattice parameters. 
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prec. 


block size p 


(N ,N h N 2 ) 


T 


MD time 


p 

1 acc 


DD 


8 4 


(4,8,12) 


1 


3000 


0.857(8) 


mass 


0.09 


(4,8,12) 


1 


3000 


0.794(8) 



Table 3: Simulation parameters for the DDHMC and the mass-preconditioned HMC algorithms. P acc denote 
the acceptance rate. 

3.4 Comparison of DDHMC and mass-preconditioned HMC 



As discussed in Sec. p.l| , it is essential for the efficiency of the DDHMC algorithm to incorpo- 



rate the multiple time scale integration scheme. It is well known that this scheme is also applicable 



to the mass-preconditioned HMC algorithm [gOj, [21[ |. We have made a direct comparison of the 
two algorithms in Nf = 2 QCD on a 16 3 x 32 lattice employing the 0(a)-improved Wilson quark 
action with the nonperturbative improvement coefficient c$w = 2.0171[f22|] and the plaquette gauge 
action at j3 = 5.2. The lattice spacing is 0.1 fm and the physical pseudoscalar meson mass is about 
600 MeV at fc ut j = 0.1355. For the mass-preconditioned HMC algorithm we employ two set of the 
pseudofermion fields which decompose the fermion determinant as 

det Q 2 = det(W" r W) det (S^) , (3.4) 

where Q is the hermitian Wilson-Dirac operator and the preconditioning operator is given by 
W = Q + p. For convenience we refer to det(W^W) as the UV part and the det(<2 2 / (W ' W)) as 
the IR part in an analogy with the DDHMC algorithm. The step sizes are chosen with the three 
integers -/Vo,i,2 in exactly the same way as the DDHMC algorithm. Simulation parameters are 
summarized in Table ||[ The block size for the DDHMC algorithm and the p parameter for the 
mass-preconditioned HMC algorithm are chosen such that H-Fcy^H are roughly the same between 
these two algorithms. This condition yields comparable acceptance ratios with M),i,2 in common. 
We employ the BiCGStab algorithm for the quark matrix inversion in both the UV and the IR parts. 

In Fig. || we plot the normalized autocorrelation function for the plaquette as a function of MD 
time. The results of both algorithms show quite similar behaviors and the integrated autocorrelation 
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Figure 4: Normalized autocorrelation functions for the plaquette (left) and the number of multiplications of 
the Wilson-Dirac quark matrix on the full lattice (right). Black (red) symbols denote the DDHMC (mass- 
preconditioned HMC) case. 
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prec. 


Tmt[P] 


T int [#mult] 


#mult 


cost[T] 


cost[#mult] 


DD 


27(10) 


22(7) 


45530(280) 


1.2(5) 


1.0(3) 


mass 


28(12) 


35(16) 


67160(380) 


1.9(8) 


2.4(1.0) 



Table 4: Integrated autocorrelation time and cost estimate for the DDHMC and the mass-preconditioned 
HMC algorithms. 



time Ti n t[P] in Table ^| are consistent within the errors. Figure |5| shows the MD-time history of the 
number of multiplications of the Wilson-Dirac quark matrix on the full lattice. The total number of 
multiplications is the sum of those required to calculate the UV and the IR forces and the Hamil- 
tonian. Their contributions are denoted by black, red, green, blue lines in order. Comparing the 
results of the DDHMC and the mass-preconditioned HMC, we observe a clear difference in the UV 
part contribution: the mass-preconditioned HMC needs more than twice of the multiplication num- 
ber for the DDHMC. This ends up in a 50% difference in the total number of multiplications. In 
Fig. ^| we also plot the normalized autocorrelation function for the total number of multiplications. 
Although the DDHMC result seems to show a slightly steeper fall-off, both results are consistent 
within the error bars. This is confirmed by the integrated autocorrelation time T; nt [#mult] in Table |j. 

Now let us compare the efficiencies of both algorithms. We define the machine-independent 
cost formula by 



cost[^] =#mult (total) /MD time x T int [^]/10 6 



(3.5) 



where the observable G is the plaquette or the total number of multiplications. In Table f| we sum- 
marize the results of cost[<^]. For both observables the DDHMC algorithm shows better efficiency 
than the mass-preconditioned HMC algorithm albeit the errors are rather large. 

There remains a couple of concerns in this study. The first one is the quark mass depen- 
dence, because our results are obtained at only one hopping parameter. The second one is the 
optimization. While we choose 8 4 block size for the DDHMC algorithm and p = 0.09 for the 
mass-preconditioned HMC algorithm since ||.Fo,i,2|| are roughly the same, these parameters may 
not be the optimal values for each of the algorithms. We leave these issues to future studies. 
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Figure 5: History for number of multiplications of the Wilson-Dirac quark matrix on the full lattice for the 
DDHMC algorithm (left) and the mass-preconditioned HMC algorithm (right). 
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Figure 6: Effective masses for the mesons (top) and the baryons (bottom) at fc u( j = 0.13727 (left) and 
0.13781 (right). Horizontal lines represent the fitting results with an error band. 
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Figure 7: Binsize dependence of magnitude of error for mesons (left) and baryons (right) at (fc u d, K s 
(0.13727,0.13640). 



4. Physical results 

4.1 Measurement of hadron masses, quark masses, decay constants 

We measure both the meson and the baryon correlators at every 10 trajectories at the unitary 
points where the valence quark masses are equal to the sea quark masses. Light hadron masses are 
extracted from a single exponential % 2 fit to the correlators with an exponentially smeared source 
and a local sink. Figure |6| shows the hadron effective masses at fc u d = 0.13727 and 0.13781 as 
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representative cases. We observe clear plateau for the mesons except for the p meson and also 
good signal for the baryons thanks to a large volume. Especially, the Q. baryon has a stable signal, 
which we use as a physical input to determine the cutoff scale later. The horizontal lines denote the 
fitting results with an error band of one standard deviation. Their widths represent the fitting ranges. 
Statistical errors are estimated by the jackknife method. In Fig. [7] we plot binsize dependence of 
magnitude of error for the mesons and the baryons at (fc u d,fc s ) = (0.13727,0.13640). For the 
pseudoscalar mesons we observe that the magnitude of error gradually increases as the bin size is 
enlarged up to about 40 MD time, beyond which it stabilizes. For other hadrons we do not find any 
clear binsize dependence. The data at other hopping parameters show similar behaviors. Based on 
this observation we choose 50 molecular dynamics time for the jackknife analysis at all the hopping 
parameters. 

We extract the bare quark mass through the axial vector Ward-Takahashi identity (AWI) by 

am « = fe 2(P(t)P(0)) ( } 

with P the pseudoscalar density and A™ p the nonperturbatively 0(a) -improved axial vector current[| 
The renormalized quark mass and the pseudoscalar meson decay constant in the continuum MS 
scheme are defined as follows: 

m m = ZA ^ +bA ^) mT , (4.2) 



f m AWI \ C 2C l 
frs = 2ku Z a \+b A ^- — P -- (4.3) 
y uq J C s p y m PS 

Here C AP are the amplitudes extracted from the correlators (A™ P (7)P(0)) and (P(t)P(O)) with an 
exponentially smeared source and a local sink, while C P is from (P(t)P(O)) with a local source and 
a local sink. The renormalization factors Zap and the improvement coefficients b A j> are evaluated 



perturbatively up to one-loop level] 24, 25] with the tadpole improvement. 



4.2 Comparison with the previous CP-PACS/JLQCD results 





lattice size 


am n 


anip 


am^i 


PACS-CS 


32 3 x 64 


0.3220(6) 


0.506(2) 


0.726(3) 


[/mini ^max] 




[13,30] 


[10,20] 


[10,20] 


CP-PACS/JLQCD 


20 3 x 40 


0.3218(8) 


0.516(3) 


0.733(4) 


[^rnin; ^max] 




[8,20] 


[9,15] 


[11,17] 



Table 5: PACS-CS and CP-PACS/JLQCD results for %, p and nucleon masses at (K ud ,K s ) = 
(0.13700,0.13640). [f m in,W] denotes the fitting range. 

We first compare the PACS-CS results on 32 3 x 64 with the previous CP-PACS/JLQCD results 
on 20 3 x 40[|I], § at (fc ud , fc s ) = (0.13700,0.13640). In Fig. § we plot the effective masses for the n 
and the p mesons. The PACS-CS and the CP-PACS/JLQCD results are consistent for the % meson, 
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Figure 8: Effective masses for the n (left) and the p (right) at (K ud ,K s ) = (0.13700,0.13640). Black and 
red symbols denote the PACS-CS and the CP-PACS/JLQCD results, respectively. 



while a slight deviation is observed for the p meson. This is numerically confirmed by the fitting 
results listed in Table ||, where we employ a single exponential % 2 fit. The nucleon mass is also 
given in Table |5[ We find 1—2% deviation for the p meson and nucleon masses, which could be 
due to possible finite size effects. 

Figure |9] shows the up-down quark mass dependence of (am n ) 2 and am p with fc s fixed at 
0.13640. For the pion mass we observe that the PACS-CS and the CP-PACS/JLQCD results are 
smoothly connected as a function of 1/ K" U( j. On the other hand, the quark mass dependence is not 
so smooth for the p meson. Although this may be attributed to finite size effects, further studies 
are needed in the p channel. 

4.3 Chiral analysis on pseudoscalar meson masses and decay constants 

We examine the chiral behaviors of the pseudoscalar meson masses and decay constants in 
comparison with the prediction of chiral perturbation theory (ChPT). Our interest exist in the fol- 
lowing points: (i) signals for chiral logarithms, (ii) determination of low energy constants in the 
chiral lagrangian, (iii) determination of the physical point with the ChPT fit, (iv) estimate of the 
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Figure 9: (am n ) 2 (left) and am p (right) as a function of l/s: ut j. Red and black symbols denote the PACS-CS 
and the CP-PACS/JLQCD results, respectively. 
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magnitude of finite size effects based on one-loop calculations of ChPT. 

We first recall the one-loop expressions of ChPT for the pseudoscalar meson masses and the 
decay constants [26] 

„,2 



1. 



2mB 1 1 + n„ - + ^ (16m(2L 8 - L 5 ) + 16 (2m + m s ) (2L 6 - L 4 )) | , (4.4) 
| = (m + m0Bo|l + ^ + ^(8('" + '«s)(2L 8 -L 5 ) + 16(2m + m s )(2L 6 -L 4 ))J(f5) 



fn = fo { 1 - 2/i s " ^ + |l (8mL 5 + 8(2m + m s )L 4 ) |> , (4.6) 



(3332} ^ 
/r = /o 1 1 - ^ - jMc " 4^ + ( 4 (™ + m s)L 5 + 8(2m + m s )L 4 ) j , (4.7) 

where m = (m u + m s )/2 and £4,5,6,8 are the low energy constants, and jltps is the chiral logarithm 
defined by 



1 m PS ,„ ( m ls 



with pL the renormalization scale. There are six unknown low energy constants Bo,/o;^4,5,6,8 in the 
expressions above. The low energy constants are scale-dependent so as to cancel that of the chiral 
logarithm (O). We determine these parameters by making a simultaneous fit for m\, m\, f n and 
h- 

We also consider the contributions of the finite size effects based on ChPT. At the one-loop 
level the finite size effects defined by Rx = (X(L) — X(oo))/X(oo) for X = m n ,mK,f-!t,fK are given 
by [0]: 

*», = ~kgi(A*)-j^gi(^), (4.9) 

R mK = U n 8i(^), (4.10) 
6 

R h = -^gi(^) - i&gi(A z ), (4.11) 

Rf K = ~^gi(^) - \$Kgx&K) - l^Sii^i) (4-12) 



with 



m ps 5 __ _ „„__ r f j\ _ £ 4m (") 



=> ps = 7IwM2 ' ^ PS = mpsL ' = L 7= K i(Vnx), ( 4 -13) 



where T^i is the Bessel function of the second kind and m(n) denotes the multiplicities in the 
expression of n = rv\ + rvl + r?. With the use of these formulae we estimate the possible finite size 
effects in our results. 

Before presenting our fitting results, it is instructive to compare the PACS-CS and the CP- 
PACS/JLQCD results for (am n ) 2 / (am^ 1 ) and f K /f n . In Fig. we plot them as a function of 

f % is normalized as 92.4 MeV in these expressions, while our results are presented in the f n = 130.7 MeV normal- 
ization. 
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Figure 10: Comparison of the PACS-CS (red) and the CP-PACS/JLQCD (black) results for {am n f / (am^ 1 ) 
(left) and fic/fx (right) as a function of am^ 1 . fc s is fixed at 0.13640. Vertical lines denote the physical 
point and star symbol represents the experimental value. 



am** 1 with K s fixed at 0.13640. The PACS-CS and the CP-PACS/JLQCD results are denoted by 
the red and the black symbols, respectively. The two sets of data together show a smooth behavior 
as a function of am™ 1 , and at K" U( j = 0.13700 (am™ 1 = 0.028) they show good consistency. It is 
important to observe that an almost linear quark mass dependence of the CP-PACS/JLQCD results 
for heavier up-down quark masses changes into a convex behavior, both for (am n ) / (am™ 1 ) and 
fx/ 'fn, as the quark mass is lowered in the PACS-CS runs. This is a characteristic feature expected 
from the ChPT prediction in the small quark mass region due to the chiral logarithm. This curvature 
drives up the ratio fx/ fx toward the experimental value as the physical point is approached. 

Let us apply the ChPT formulae (gg)-(g^) to our results at four points (jQ d , jc s ) = (0.13781,0.13640), 
(0.13770,0.13640), (0.13754,0.13640), (0.13754,0.13660). For these points, the p meson mass sat- 
isfies the condition that m p > 2m n . The measured bare AWI quark masses are used for m and m s in 
eqs.([0|)-(|4/7|). The heaviest pion mass at (fc ud , K s ) = (0.13754,0.13640) is about 430 MeV with 
the use of the cutoff determined below. We summarize the pion masses and the unrenormalized 
AWI quark masses in Table ^. The fit results are shown in Fig. 1 1 , where the black solid lines 
are drawn with fc s fixed at 0.13640 and the black dotted lines are for K s = 0.13660. The red solid 
symbols represent the extrapolated values at the physical point whose determination is explained in 



Kud 




am n 




am™ 1 


0.13700 


0.13640 


0.32196(62) 


0.02800(20) 


0.04295(30) 


0.13727 


0.13640 


0.26190(66) 


0.01895(13) 


0.04061(18) 


0.13754 


0.13640 


0.18998(56) 


0.01020(11) 


0.03876(18) 




0.13660 


0.17934(78) 


0.00908(7) 


0.03257(17) 


0.13770 


0.13640 


0.13591(88) 


0.00521(9) 


0.03767(10) 


0.13781 


0.13640 


0.08989(291) 


0.00227(16) 


0.03716(20) 



Table 6: Pion masses and unrenormalized AWI quark masses. 
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Figure 11: Fitting results for (ara,) 2 /(flffi™) (left) and fx/ fit (right)- Red solid (open) symbols denote 
the extrapolated values at the physical point by the ChPT formulae without (with) the finite size effects. 



Sec, below. The heaviest point at (rC ud ,rC s ) = (0.13754,0.13640) is not well described by ChPT 
both for (am K ) 2 / \am^ 1 ) and fx/ fn, and % 2 l&.o.f. is rather large (see Table [7]). 



L,(jU =m n ) 


PACS-CS 


PACS-CS with FSE 


exp. value[28 


] MILC[ 


29] 


U 

u 

2L(, — L4 
2U-L5 


0.25(11) 
2.28(13) 
0.16(4) 
-0.59(5) 


0.23(12) 
2.29(14) 
0.16(4) 
-0.60(5) 


0.27 ± 0.8 
2.28 ±0.1 

0± 1.0 
0.18 ±0.5 


0.1(2)(2) 
2.0(3)(2) 
0.5(1)(2) 
-0.1(1)(1) 


X 2 /d.o.f. 


2.1(1.4) 


2.1(1.4) 









Table 7: Results for the low energy constants together with the phenomenological estimates 
MILC results 



and the 



The results for the low energy constants are presented in Table [7| where the phenomenological 
values with the experimental inputs[^|] and the MILC results^] are also given for comparison. 
The renormalization scale is chosen to be = 0.547 GeV. For L4 and L5 governing the behavior of 
fn, fx, our results show good agreement with both the phenomenological estimates and the MILC 
results. On the other hand, some discrepancies are observed between three results for 2L& — L4 and 
2L% — L5 which enter into the ChPT formulae for m\ and m 2 K . 



In Fig. |11| we also draw the ChPT fit results including the finite size effects. The green solid 
lines are drawn for fc s = 0.13640 and the green dotted ones for K s = 0.13660. The fit curves 
with and without the finite size effects are almost degenerate for am^ 1 > 0.003, but deviations 
appear closer to the physical point, for which the extrapolated values are plotted by the open and 
solid red symbols. This feature is understood by Fig. [j/J where we plot the magnitude of Rx for 
X = m K ,mK,fit,fK with L = 2.8 fm as a function of m % ( we note that R mps > and Rf ps < 0). The 
finite size effects are less than 2% for mps and fp$ at our simulation points. For mps this is true 
even at the physical point, while for f n the finite size effects cause the value to decrease by 4%. 
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Figure 12: \Rx\ (R mps > and Rf ps < 0) forX — m K ,mK,fx,fK with L = 2.8 fm as a function of m n . Solid 
vertical line denotes the physical point and the dotted ones are for our simulation points. 



4.4 Physical point and light hadron spectrum 

In order to determine the up-down and the strange quark masses and the lattice cutoff we need 
three physical inputs. We try the following two cases: m K ,m K ,mci and m 7C ,m K ,m ( j ) . The choice of 
mci has theoretical and practical advantages: the Q. baryon is stable in the strong interactions and its 
mass, being composed of three strange quarks, is determined with good precision with small finite 
size effects. We also choose for comparison. We employ the NLO ChPT formulae for the chiral 
extrapolations of m K , niK, f% and fy. A simple linear formula nihad = Co + c\ • '"ucT 1 + c 2 ' m ^ Nl 
is used for the other hadron masses, employing data in the same range K u a > 0.13754 as for the 



pseudoscalar mesons. In Fig. 13 we show the linear chiral extrapolations for ma and ma- The solid 
lines are drawn with fc s fixed at 0.13640 and the dotted ones are for fc s = 0.13660. We observe 
that the quark mass dependences for m<f, and mq at K u a > 0.13754 are well described by the linear 
function. 

The results for the quark masses and the lattice cutoff are listed in Table |8[ where the errors are 
statistical. The two sets of results are consistent within the error. The quark masses are smaller than 
the recent estimates in the literature. We note, however, that we employed the perturbative renor- 
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Figure 13: Linear chiral extrapolation for aran (left) and am^ (right). Solid (dotted) lines are drawn with 
K s — 0. 1 3640 (0. 1 3660). Red open symbols denote the extrapolated values at the physical point with a linear 
form. 
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input 


a~ l [GeV] 


m^ s [MeV] 


mf s [MeV] 


A 


h 


fa/fx 


ma 
m$ 


2.256(81) 
2.248(76) 


2.37(11) 
2.38(11) 


69.1(25) 
69.4(25) 


144(6) 
143(6) 


175(6) 
175(5) 


1.219(22) 
1.219(21) 



Table 8: Cutoff, renormalized quark masses, pseudoscalar meson decay constants determined with mjj and 
m^ inputs. 
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Figure 14: Light hadron spectrum extrapolated at the physical point with £2-input (red) and 0-input (blue). 
Horizontal bars denote the experimental values. 



malization factors to one-loop level which may contain a sizable uncertainty. A nonperutrbative 
calculation of the renormalization factor is in progress using the Schrodinger functional scheme. 

In Table || we also present predictions for the pseudoscalar meson decay constants at the phys- 
ical point using the physical quark masses and the cutoff determined above, which should be com- 
pared with the experimental values f % = 130.7 MeV, f K = 159.8 MeV, f K /f K = 1.223. A 10% 
discrepancy in the magnitude of f n and fg might be due to use of one-loop perturbative Za since 
the ratio shows a good agreement. A nonperturbative calculation of Za is also in progress. 

In Fig. [14| we compare the light hadron spectrum extrapolated to the physical point with the 
experiment. The results for the Q-input and the <p -input are consistent with each other, and both 
are in agreement with the experiment albeit errors are still not small for some of the hadrons. This 
is an encouraging result. However, further work is needed since cutoff errors of 0((ciAqcd) 2 ) are 
present in our results. 



5. p-nn mixing 

Since our simulations are carried out at sufficiently small up-down quark masses, it would be 
interesting to investigate the p-%% mixing effects. We find that the rest mass m p is always smaller 
than the two-pion energy 1\Jm\ + (2n/L) 2 for all the hopping parameters, and hence the p meson 



at rest cannot decay into two pions. However, as illustrated in Fig. |15j, for a moving p with a unit of 



momentum, i.e.,, its energy yw?p + (2n/L) 2 becomes larger than the energy of a moving pion and 

a pion at rest given by \jm\ + {2%/L) 2 +m n when the up-down quark mass is sufficiently reduced. 

Let us consider two types of the p meson propagator with the momentum 2n/L: pn(2n/L) 
with polarization parallel to the spatial momentum and p±(2n/L) with polarization perpendicular 
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Figure 15: Energy levels of the p meson and the two pion states without the total momentum p = (left) 
and with p = p m \ n = 2n/L (right) as a function of the up-down quark mass. K s is fixed at 0. 13640. 




to the spatial momentum. Phenomenologically the p-%% coupling is described by g pK 7i^abcP C j!i^ b d^7i c , 
which favors p|| (2k /L) ->■ k(2k / L)k(<S) to p L (2n/L) -> k(2k / L)k(0) . We expect that the p\\ (2%/L) 
propagator is more strongly affected by the mixing effects than the p±(2n/L) correlator. Since the 
mixing effects push up the upper energy level further and push down the lower energy level as 
shown in Fig. [if], they could be detected by measuring the R function defined by 

(P\\(P,t)p}(P,0)) large t -(E n -E„ )t 
R(A = — 1 L, Ze 1 "II . (5.1) 

(p ± (p,t)pl(p,0)) 



In Fig. 16 we plot log | R(t)\ as a function of t. The dotted horizontal lines denote R(t) = (E/m p ) 2 
, which is determined kinematically in the mixing-free case. The solid lines represent the fitting 
results with a single exponential form over 5 < t < 1 1 . The data show clear positive slopes which 
indicate E p ,, < E p± . We also observe that the magnitude of the energy difference is rather small for 
K"ud < 0.13754, while it grows rapidly as the up-down quark mass is reduced for fc u d > 0.13754. 
This feature may suggest that the (p|| (p,t)p^(p,0)) correlator is getting dominated by the %% state 
toward the smaller up-down quark masses. In order to obtain a definite conclusion, we need more 
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detailed investigations with increased statistics. 
6. Summary 

We have presented a status report of the PACS-CS project which aims at a 2+1 flavor lattice 
QCD simulation toward the physical point. With the aid of the DDHMC algorithm for the up- 
down quarks we have reached m n = 210 MeV, which roughly corresponds to m^(jx = 2 GeV) = 
5.6 MeV, on a 32 3 x 64 lattice using the 0(a) -improved Wilson quarks. Thanks to the enlarged 
volume compared to the previous CP-PACS/JLQCD work, we obtain good signals not only for the 
meson masses but also for the baryon masses. Our results for the hadron spectrum at the physical 
point show a good agreement with the experimental values. 

At present we have just started the simulation at the physical point. We are also calculating the 
nonperturbative renormalization factors for the quark masses and the pseudoscalar meson decay 
constants in order to remove perturbative uncertainties in these important quantities. Once these 
calculations are accomplished, the next step is to investigate the finite size effects at the physi- 
cal point, and then to reduce the discretization errors by carrying out calculations at finer lattice 
spacings. 

Acknowledgment 

We would like to thank all the collaboration members for discussions and, in particular, 
A. Ukawa for a careful reading of this report. Numerical calculations for the present work have 
been carried out under the "Interdisciplinary Computational Science Program" in Center for Com- 
putational Sciences, University of Tsukuba. This work is supported in part by Grants-in-Aid 
for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology 
(Nos. 13135204, 15540251, 17340066, 17540259, 18104005, 18540250, 18740139). 

References 

[1] CP-PACS and JLQCD Collaborations, T. Ishikawa et ai, PoS (LAT2006) 181. 

[2] CP-PACS and JLQCD Collaborations, T. Ishikawa et ai, hep-lat/0704. 193. 

[3] CP-PACS and JLQCD Collaborations, S. Aoki et ai, Phys. Rev. D73 (2006) 034501. 

[4] Y. Iwasaki, preprint, UTHEP-1 18 (Dec. 1983), unpublished. 

[5] JLQCD Collaborations, S. Aoki et ai, Phys. Rev. D65 (2002) 094507. 

[6] PACS-CS Collaboration, S. Aoki et ai, PoS (LAT2005) 111. 

[7] PACS-CS Collaboration, A. Ukawa et ai, PoS (LAT2006) 039. 

[8] PACS-CS Collaboration, Y. Kuramashi et al. PoS (LAT2006) 029. 

[9] M. Luscher, JHEP 0305 052 (2003); Comput. Phys. Commun. 165 (2005) 199. 

[10] A. Kennedy, Nucl. Phys. B (Proc. Suppl.) 140 (2005) 190. 

[11] PACS-CS Collaboration, K-I. Ishikawa et al. PoS (LAT2006) 027. 

[12] PACS-CS Collaboration, N. Ukita et ai, PoS (LAT2007) 138. 



19 



Nf = 2+1 dynamical Wilson quark simulation toward the physical point 



[13] PACS-CS Collaboration, D. Kadoh et al, PoS (LAT2007) 109. 

[14] L. Del Debbio et al, JHEP 0602 (2006) 011. 

[15] L. Del Debbio et al, JHEP 0702 056 (2007); JHEP 0702 (2007) 082. 

[16] J. C. Sexton and D. H. Weingarten, Nucl. Phys. B380 (1992) 665. 

[17] M. Luscher, Comput. Phys. Commun. 156 (2004) 209. 

[18] A. Ukawa, Nucl. Phys. B (Proc. Suppl.) 106 (2002) 195. 

[19] CP-PACS Collaboration, A. Ali Khan et al, Phys. Rev. D65 (2002) 054505. 

[20] M. Hasenbusch, Phys. Lett. B519 (2001) 177. 

[21] M. Hasenbusch and K. Jansen, Nucl. Phys. B659 (2003) 299. 

[22] ALPHA Collaboration, K. Jansen and R. Sommer, Nucl. Phys. B530 (1998) 185. 

[23] CP-PACS/JLQCD and ALPHA Collaborations, T. Kaneko et al, JHEP 0704 (2007) 092. 

[24] S. Aoki et al, Phys. Rev. D58 (1998) 074505. 

[25] Y. Taniguchi and A. Ukawa, Phys. Rev. D58 (1998) 1 14503. 

[26] J. Gasser and H. Leutwyler, Ann of Phys. 158 (1984) 142; Nucl. Phys. B250 (1985) 465. 

[27] G. Colangelo, S. Diirr and C. Haefeli, Nucl. Phys. B721 (2005) 136. 

[28] G. Amoros, J. Bijnens and P. Talavera, Nucl. Phys. B602 (2001) 87. 

[29] C. Bernard et al, arXiv:hep-lat/061 1024. 



20 



